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ABSTRACT 

This paper discusses numerical solutions of a hyperbolic initial boundary value 
problem that arises from acoustic wave propagation in the atmosphere. Field equa- 
tions are derived from the atmospheric fluid flow governed by the Euler equations. 
The resulting original problem is nonlinear. A first order linearized version of the 
problem is used for computational purposes. The main difficulty in the problem as 
with any open boundary problem is in obtaining stable boundary conditions. Approxi- 
mate boundary conditions are derived and shown to be stable. Numerical results are 
presented to verify the effectiveness of these boundary conditions. 
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1. Introduction 

In this paper we present results concerning a model problem that arises in the study of nonlinear 
acoustic wave propagation in the atmosphere. The derivation follows after Cole & Greifinger [1], and 
the details are presented in the appendix. In an earlier study Hariharan [3,4] discussed a strictly two 
dimensional version of this problem. In these references the well-posedness of the problem was dis- 
cussed, and results concerning stable boundary conditions were presented. The present paper discusses 
an axisymmetric three dimensional model. It is well known that the nature of wave propagation varies 
drastically between two and three dimensions and two dimensional results are not special cases of the 
three dimensional results. To this end the present study is motivated by possible comparisons with 
future experimental results. Needless to say the experiments will be three dimensional, and any com- 
parable theoretical model also has to be in three dimensions. Most of the acoustic wave propagation 
experiments are done in a laboratory environment; therefore, simulation of actual nonlinear effects is 
rather difficult. In particular, long range wave propagation in the atmosphere is a difficult situation to 
simulate in the laboratory. As a result we shall confine ourselves to obtaining numerical results for only 
a linear version of the problem, which is a difficult problem in itself. Moreover, a careful study of this 
problem will lead to meaningful nonlinear corrections later once we have some experimental results. 
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In [2], it was pointed out that the nonlinear problem can be solved by obtaining a sequence of 
linear problems successively. Each succeeding problem contains the solution of the previous problem 
as forcing terms. It was also proved in [2] that the problems were well-posed. We develop analogous 
arguments for the axisymmetric three dimensional case in this paper. However, our main focus here is 
to obtain numerical solutions for the first linear problem that arises in the present study. 

Let us present the problem (see appendix) under consideration. The model assumes the following 
situation as explained in the appendix. 

(i) Ground (earth) is flat. 

(ii) Atmosphere is isothermal. 

(iii) Pressure and density vary exponentially (standard model). 


z 



z 

b 


| source 


Figure (1.1 b) 




The problem is posed on a half space (-<*> ,oo)x(-oo ,«>)x [0,<>o ) with a point source (figure 1.1 a) at z 0 
emitting acoustic waves. Due to axisymmetry, we consider the problem in the quarter space 
[0,°° )x [0,°o) (in polar coordinates). If p,u,v, and p denote the acoustic density, acoustic velocity com- 
ponents in the r and z directions and acoustic pressure respectively, then the governing equations are: 


( 1 . 1 ) 
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p, + u r + — + w z - w = 0 


u t + - Pr = o 
w t + — p z - = 0 


Pt + YUr + yw z + Yj - W = f(r,z,t) 

where y = 1.4, f(r,z,t) describes the nature of the acoustic source and r = V x 2 +y 2 . 
Boundary Conditions are: 


( 1 . 2 ) 

(1.3) 

(1.4) 


Initial Conditions are: 


w(r,0,t) = 0, t> 0 
u(r,z,t)/r is bounded as r— >0. 


(1.5) 

( 1 . 6 ) 


p = p = u = w = 0 for t = 0. (1.7) 

For convenience, we shall call the initial boundary value problem defined through (1.1) - (1.7) 
problem (P). As mentioned earlier, problem (P) is the first linearized problem that results from the 
acoustic equations. Our goal is to solve (P) numerically. As in any open boundary problems, the prob- 
lem at hand cannot be solved numerically unless the infinite region is truncated to a finite one. When 
we do such a truncation, we also have to provide boundary conditions on these truncated boundaries. 
They have to simulate the behavior at infinity and should be absorbing boundary conditions. Moreover 
the problem in the truncated region must be well-posed. Thus the plan of the paper is as follows: First 
we construct radiation boundary conditions that simulate the behavior at infinity at finite distances. 
Then we show that the problem in the truncated region is well-posed. Then we present the numerical 
technique along with some numerical boundary conditions that are required to stabilize the numerical 
scheme. Finally, we present some examples to validate the numerical scheme and to interpret the phy- 
sics behind the problem. 


2. Absorbing boundary conditions 

As mentioned earlier, it is essential to truncate the boundaries that are at infinity to the ones at 
finite distances. Thus we do the truncation into a finite region [0,L] x [0,H] s Q. The situation is 
shown in figure 2.1. 
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Figure 2.1 ^ ^ 

To obtain radiative or absorbing boundary conditions, there are a variety of procedures. One way to 
construct such boundary conditions is to obtain appropriate inflow and outflow variables for the system. 
This procedure is similar to obtaining appropriate Riemann invariants. For successful implementations 
in one dimension see Hariharan and Lester[4]. However, the problem under consideration is a quasi 
two-dimensional problem. An exact construction of the Riemann invariant is not possible. But an 
acceptable procedure which is approximate is to consider one-way wave equations in the r and z direc- 
tions and obtain the Riemann invariants. To do this let us note that the system of equations under con- 
sideration (equations (1.1) -(1.4)) can be written in the form: 

u t + Au r + Bu z + Cu = F (2.1) 

where, 

0 10 0 
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F = [0,0,0,f] T and u = [p,u,w,p] T - 

At large distances the effect of the source term vanishes. Thus the one wav wave equation in the r 







direction is 
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u t + Au r = 0 (2.2) 

and in the z direction it is 

u t + Bu z = 0. (2.3) 

We note that we have neglected the lower order terms in this process, and this again is an acceptable 

procedure to obtain boundary conditions. The final step is to diagonalize the matrices A and B through 

a change of variable u of the form u = Tv, so that T _1 AT and T _1 BT are diagonal. The elements of 

the resulting diagonal matrices will contain the eigenvalues of A and B. The matrix T in each case is 

constructed using the eigenvectors of A and B. This leads to the construction of the inflow and outflow 

variables at the boundaries Tj and r 2 (figure 2.1). The inflow variables on these boundaries are p - yu 

on Tj and p - yw on T 2 . Thus for radiating boundaries we require that these variables be zero and 

obtain the radiation boundary conditions. They are 

p - yu = 0 on n (2.4) 

and 

p - yw = 0 on r 2 . (2.5) 

In the next section we shall show that the problem (P) together with boundary conditions (2.4) and 

(2.5) yield a well-posed problem. 

One should note that these boundary conditions are the lowest order boundary conditions, and as 
such one has to keep the truncated boundaries T x and T 2 at far distances to obtain minimal reflections 
from infinity. To obtain higher order boundary conditions, one has to resort to a different strategy. One 
needs to obtain the dispersion relations of the problem under consideration. These dispersion relations 
can also be interpreted in the pseudo differential operator terminology. Whatever the case may be, the 
ideas lead to obtaining boundary operators that dictate radiativity. For associated details we refer readers 
to [5,6]. In these references only simple wave equations or problems that are not of the characteristic 
type are discussed. The current problem as we will see has a complicated structure. To investigate 
further, let us consider the governing equations of (P). Suppose we are interested in obtaining radia- 
tion conditions on the boundary Tj . We take the Laplace transform in time of equations ( 1. 1) -( 1 .4) to 


- 6 - 


obtain: 


_ _ u _ 

sp + u r + — + w z - w= 0 

- 1 - „ 
su + — p. = 0 

y 


sw + — p z - -5* — = 0 

y y 

sp + yu t + ywj + y~ - w = 0. 


Here for any function h(r,z,t) the Laplace transform is defined by 


which will be inverted according to 


h(r,z,s) = J 0 e st h(r,z,t) dt 

) 

h(r,z,t) = — ~ f . e st h(r,z,s) ds. 
2n i J ~ lo ° 


( 2 . 10 ) 


( 2 . 11 ) 


( 2 . 12 ) 


Since our goal is to obtain radiating boundary conditions on r 2 , i.e., in the z direction, we will 
transform the r variable also to obtain a system of ordinary differential equations in the z direction. Due 
to the presence of cylindrical symmetry, it is natural to use the Hankel transforms. For any function 
g(r), the Hankel transform and its inverse are defined respectively by 


g(w) = J Q g(r)J y (cor)r dr ( y > - j) 

oo 1 

g( f ) = J 0 g(co)J y (cor)co dto (y > 


(2.13) 


(2.14) 


(p,w,p)(z;w,s) = J Q (p,w,p)(z,r,s) J 0 (tor)r dr 


(2.15) 


u(z;co,s) = J o u(z,r,s) Jj(o)r)r dr. 


(2.16) 


With these definitions we apply the Hankel transforms to equations (2.7)-(2.10). We note that the 
definition of the transform is different for u than the other dependent variable. This is due to the singu- 
larity in r. However, the recurrence relations of the Bessel functions J 0 and Jj account for the singular- 
ity and the final transformations look as if we took the Fourier transform ignoring the singular terms. 
We shall show the computations to one equation. 
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Multiplying equation (2.7) by J 0 (wr)r and integrating over r on the interval [0,«), we get 

sp + J Q (u r + Jo(cor)r dr + w z - w = 0. 

We rearrange the integral term and integrate by parts to get 


J 0 (ru)rJo(wr) dr=-J Q u[J 0 ((or)]'r dr. 
The recurrence relation needed here is 


(2.17) 


[Jo((or)]'= -concur). 

Thus the right hand side of (2.17) becomes cou. Similar calculations applied to the remaining equations 
yield the following system of ordinary differential equations in z: 


sp + cou + w z - w = 0 

(2.18) 

su + — f>= 0 
Y 

(2.19) 

sw + — p z - = 0 

(2.20) 

Y Y 


sp + you + yw z - w = 0. 

(2.21) 


Seeking solutions of the form e Xz , one obtains the following characteristic relations for the above sys- 
tem: 


where, 


and where, 


X = 



H(s,co), 


lt(s,to) = 


-7 + S 2 + CO 2 + 
4 


(3 2 co 2 /s 2 


( 2 . 22 ) 


(2.23) 


p 2 =(Y-l)/y 2 (> 0). 

This relation is the same as the dispersion relation for the strictly two dimensional model. It was argued 
in [2] that for propagation the ratio — must be smaller than one, and the outgoing waves are given by 
(2.22) with the sign. These lead to the observation that the dispersion relation 

^ = \ - M-(s.co) (2.24) 


can be approximated for large s and the approximations arising from (2.23) are found in [2], The first 



two approximations are 


- 8 - 



To obtain the absorbing boundary operators, we multiply (2.25) by e Xz 
variable, in this case the acoustic pressure. This leads to the differential 
transform domain 


(2.25) 

(2.26) 

and an appropriate dependent 
operator (approximate) in the 


Pz={e— sf>. (2.27) 

Now we apply the inversion formulas (2.12) and (2.14) to get the required approximate boundary con- 
dition 


Pz = -JP - Pt- (2.28) 

For the second approximation (2.26), to carry out this procedure one must multiply the equation by s. 
This requirement is easily seen from the equation because it contains s in the denominator, which 
causes difficulties in the inversion. Then the rest of the procedure holds, and the resulting boundary 
operator is 

Pzt = yPt - Ptt - -jj- - jPrr. (2.29) 

A similar procedure can be used to derive boundary conditions on the boundary Tj. However, from 
(2.29) one must note that the boundary condition is second order, while the system we want to solve is 
first order. This leads to an undesirable situation for numerical implementations. Also, we have not 
proved that (2.28) or (2.29) are stable boundary conditions, i.e., they yield a well-posed problem in the 
sense of an energy estimate discussed in the next section. However, the above higher order boundary 
conditions are possible candidates to explore computational efficiency. 

It is interesting to note that in equation (2.28) if we neglect the lower order term we obtain a 
further approximation in the form: 

Pz=-p t . (2.30) 

Moreover, if we incorporate equation (L3) neglecting the lower order terms, equation (2.30) becomes 



(2.31) 
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Pt~ Y w t = 0. 

Integrating (2.31) using the initial conditions, one obtains the approximate boundary condition 

p - yw = 0 (on r 2 ) (2.32) 

which is the same as (2.4) which we obtained through the characteristics. For this reason we shall refer 
to (2.32) as the zeroth order boundary condition and (2.28) as the first order boundary condition. Our 
discussion and numerical experiments on the problem pertains only to the zeroth order condition. 
There are two reasons for such a consideration. First is the difficulty in establishing well-posedness of 
the problem with the higher order boundary conditions. The second one is due to the fact that it is 
desirable to establish a solution procedure with a simpler set of boundary conditions, which has the 
theoretical backing rather than concern ourselves about the accuracy and computational efficiency with 
higher order conditions, for which we do not have the analysis yet. 


3. Energy Estimates 

Here, we show that the problem (P) together with the absorbing boundary conditions (2.32) lead 
to a well-posed problem. To do so we consider an appropriate energy for the problem in the truncated 
region Q and show it is bounded by the energy supplied through the source. Recall that the governing 
equations can be written in the form (equation (2.1)) 

u t + Au r + Buj+ Cu = 0. (3.1) 

Also, physical considerations yield the boundary conditions on T 3 and on T 4 (see figure 2.1) 


u l r = 0 (on the axis) 


(3.2) 


and 


wlj- 4 =0 (on the ground). (3.3) 

Radiation conditions are 

(p-Yu)l ri =0 (3.4) 

(p - yw) lf 2 = 0. (3.5) 

The definition of the vectors u, F , and the matrices A, B , and C are as in section 2. Also we have the 


initial condition 
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u = 0 (t=0). (3.6) 

We call the problem defined through (3.1) - (3.6) (P0- We show here the following: 

Theorem 3.1 

Problem (P') is well-posed satisfying an energy estimate of the form 

J o I le _T,t u I Iq dt < Ctt’J,, I le _T,t F llo dt, (3.7) 

where the norm 1 1. 1 1 0 is the L 2 norm in Cl , C is a constant independent of F, for all rj > K, a positive con- 
stant . 

To prove this theorem we note that the system (3.1) is a system of hyperbolic partial differential 
equations, but it is not strictly hyperbolic. As such Kreiss’s stability theory [7] for the given initial 
boundary value problem is not directly applicable. However, we obtain energy estimates as in theorem 
3.1. To do that we need a matrix G such that the matrices M = GAG” 1 and N = GBG -1 are symmetric. 
If we multiply equation (3.1) on the left by the matrix S = GG T , then it takes the form 

Su t + Pu r + Qu z + SCu= SF (3.8) 

where P = SA and Q = SB are symmetric and S is positive definite. The matrix S is called the sym- 
metrizes 

Lemma 3.1 

There exists a matrix G such that the matrices M = GAG” 1 and N = GBG” 1 are symmetric. 

Proof of lemma 3.1 can be seen in reference^]. The basic idea of construction of G relies on diagonal- 
izing the A matrix using the matrix formed by its eigenvectors and multiplying the transformed system 
by a diagonal matrix whose entries are unknown [8]. These entries are determined so that the 
transformed B matrix is symmetric. If T is the matrix formed by the eigenvectors of A , then T _1 AT 
will be diagonal but not T^BT. Further we use a diagonal matrix D = (a,p,y,S) and form 

D” 1 T" 1 ATD and d” 1 T" 1 BTD. 

The first matrix will still be diagonal so that it is symmetric while the second one becomes symmetric 
for the choice of a= 1, P= V 2, 1 , and 8= 1. With this choice of parameters we obtain 
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s = 


0 0 -■ 


0 j 0 

0 0 \ 

-—0 0 
Y 


0 0 0 
0 0 0 


p = 


Q = 


Y 

0 

0 

3 

2y 2 

0 

1 


2y 

oooo 

o o o 

2y 


0 0 0 
0 0 0 

0 0 0 

1 


sc = 


0 

0 

1 

2y 

0 


0 0 

0 

0 

0 

1 

2yr 


0 

0 

1 

2y 

o 


2y 

- 1 + — 
Y 

0 

0 

J_ 3_ 

y 2y 2 


0 

0 

1 

2y 

o 


Returning to the proof of the theorem we now multiply (3.8) on the left by u T to obtain 

u T Su t + u T Pu r + u T Qu z + u T (SC)u = u t SF. 

Since the matrices S, P , and Q are symmetric, we can write (3.9) as 


•j(u T Su) t + -j(u T Pu) r + •~(u t Qu) z + u T SCu = u t SF. 

Let us define a quantity H(t) that measures the energy associated with the system (3.10): 

H L 

H(t) = J ^(u T Su) dx = 27tj o J Q r(u T Su) dr dz. 

Using (3.10) we obtain, 

dH T , H , L r H r L t 

“dT = 271 L /o Jo r <“ Pu >' dr dz - Jo Jo r(u Qu)! dr dz 

- 2j o j o t u T SCu dr dz + 2 J o J Q ru T SF dr dz j . 


(3.9) 


(3.10) 


(3.11) 


Applying the divergence theorem, it follows 

IT = 271 [" J r 1 LuTpu dz + J 0 H J o uTpu dr dz - J r/ uT Q u dr 

f r H L 

+ J r ru T Qu dr - 2j q J Q ru T (SC)u dr dz 

+ Jo Jo uTsFrdrd z]. 

It is easy to see that, 

- 2 Jo Jo ruT ( sc -|r> udrdz * K[j o H J o L ru T Sudrdz_ 

where 

0 0 -1+— o 

Y 

0 0 0 - — 

4yr 

7~0 0 --^-y 

ly 2 ' 

0 "T" "--A - 0 I 

4yr y 2y 2 

and K is a positive constant. 

Thus, 

dH r 

— < 2n |- J ri Lu T Pu dz - J ru T Qu dr + J p ru T Qu dr 

H L 4 
+ 4jtJ o J o uTSFr dr dz + KH. 

The boundary condition w = 0 on gives 



J r ^ru T Qu dr = 0. 

Now boundary conditions on T] and r 2 yield 

J r Lu t Pu dz 
J r ^ru T Qu dr 

Then (3.14) takes the simple form 

dH H L 

— < 4tc I I u T SFr dr dz + KH. 
dt J o J o 


= r H L i>iL dz > o 

Jo y 

= jAp — dr> 0. 
Jo y 


(3.12) 


(3.13) 


(3.14) 


(3.15) 


Moreover, 
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By the Schwartz inequality, 

|(Gu) t (GF) | < -2-u t Su+ -^F t SF - 

Hence 

ll.”jo L “ TSFr<lr ' ,I l S lJ. H J» t “ Ts “‘ lrd2+ 

for all t\ > 0. Thus we obtain 

f- S (n+K)H + ^-/ 0 "j o VsFrdrd Z 

or 

JL( e -Oi+K)t H) < ^L e -(ti+ K )tJ o H J o L F T S Fr dr dz. (3.16) 

Integrating with respect to t and using the initial condition u = 0 at t = 0 gives us 

27te _(11+K)x J 0 H ] 0 L u T Sur dr dz < ~j o ”j o L e - (T1+K)t F T SFr dr dz , (3.17) 

Finally, we integrate with respect to t to obtain the required energy estimate 

J J I le' T|, u 1 1 0 2 dt < J J I le _T,t F 1 1 0 2 dt 

where C is a constant independent of F, for all T| > K. 

4. Numerical Considerations and results 

The numerical method used here is a straightforward explicit method, second order accurate in 
space and time, commonly known as MacCormacks scheme. A split form of this scheme has been used 
for numerical calculations. Recall that the system under consideration can be written in the form 

u t + Au, + Bu z + Cu = F. (4.1) 

We further write the system in the form 

u t + F r + G z = H, (4.2) 

where H = F - Cu. We decompose H into H[ and H 2 so that equation (4.2) is split according to 

u t + F t = H t (4.3) 

u t + G z = H 2 . (4.4) 

Let us denote the solution operators of (4.3) and (4.4) by L r (At) and L z (At). Thus the solution of 
equation (4.2) is advanced according to 


u n+1 = L r L z u' 


(4.5) 
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In the calculations reported here H! and H 2 are chosen as 


H ; = H 2 = “7 H = •— w - — , 0 , P ■ P , f + w - X!L 

1 2 2 r Y r 


Corresponding flux quantities F and G are 


F = u , — p , 0 , yu 

y 


G = w , 0 , — p , yw . (4.8) 

y 

Note that quantities Hj and H 2 have singular terms near the origin. They are handled as follows: at 
r = 0, — is replaced by u, using L’Hospital’s rule. Thus at r = 0, F and H t are modified according to 


_ 3 1 A 3 

F= yu , — p , 0 , — Yu 


H = w , 0 , •£-£- , f + w . 

y 


(4.10) 


Moreover, u = 0 on the axis; thus (4.9) takes the form F = |^0 , — p , 0 , 0J on the axis. Absorbing 

boundary conditions (2.4) and (2.5) are applied on Tj and r 2 . Also characteristic extrapolations were 
applied on these boundaries. For examples, on T x the radiative boundary condition is p- yu = 0, 
which is an incoming characteristic (set equal to zero). In each cycle we compute p + yu = T, the out- 

T T 

going characteristics. Solving for p and u from these two relations, we obtain p = — and u = 
Similar consideration was given on T 2 . 

Numerical computations were performed on a 41 x 21 grid (i.e., 41 point in the r direction and 21 
in the z direction). Nondimensional lengths L and H were chosen to be 30 and 10 respectively. Higher 
resolution in the z direction was essential for the stability of the code. At was chosen sufficiently small 
to meet the stability criteria of the scheme. 

Special consideration on the choice of the source was given as indicated below. The source was 
modelled according to 


f(r,z,t) = 


cos t 
An/ Az 
0 


(r,z) = 0 


(4.11) 
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which is in l 2 while having a property that is closely related to ^ — 8 ( 2 ) cos t. Physically such a source 

corresponds to a point source which oscillates harmonically. Effectiveness of the code was verified by 
comparing a solution that is reported in Cole and Greifinger [ 1 ] for a point source both in space and 
time of the form 


f(r,z,t) = ^- 8 (z) 8 (t). 

This was modelled by an l 2 source for numerical purposes as 


f(r,z,t) = 


1 


(4.12) 


(4.13) 


ArV AzAt 

An interesting case here occurs when the source takes the form (4.11). According to our esti- 
mate given by theorem (3.1), it is valid only for finite time. As time steadily grows the solution may 
become unbounded. Our computational results indicate such a phenomenon. Up to five periods (lCbr 
units) of time, the solution propagated through the radiation boundaries smoothly. Above that time 
the solution had a tendency to grow in the computational domain. We shall show this result up to ten 
periods of time. 

All the results presented herein correspond to a source located at r = 0, z = 0 as indicated by 
(4.11). In figure (4.1) we show the time history of the pressure wave on the radiation boundary at the 
sixth grid point in the z direction. As we see from the figure, the solution starts from a state of rest 
and reaches a harmonic state. In particular the amplitude diminishes. For this case we have also plot- 
ted the time history of the wave at the same z level at four different locations of the computational 
domain: r = 5, r = 10, r = 15 and r = 20 . The radiation boundary T x is at r = 30 (see figures (4.2) 
(a) -(d)). 

The sound pressure level distribution is of interest from the physical point of view. This is illus- 
trated in figure (4.3). In the absence of gravitational effects, the sound pressure level will be monoton- 
ically decreasing since the pressure decay is inversely proportional to the distance of observation. 
Clearly this is not the case as seen in figure (4.3), This is one of the features of the propagation of 
sound in the atmosphere. At certain regions sound may not be detected while at the same time, it is 
possible to detect the sound at a further distance from the source. 
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Figures (4.4) and (4.5) repeat the calculations of the above situation for ten periods of time. Fig- 
ure (4.4) shows that roughly after seven periods of time, the solution begins to grow. Further increases 
in time lead to more growth in the solution. This occurs even with higher order corrections in radiation 
boundary conditions. Figure (4.5) illustrates the sound pressure level distribution in this situation. 
The phenomena for growth as well as nonlinear correction to the problem shall be considered in a later 
paper. 
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Appendix 


Derivation of the field equations 


We shall begin with the statement of the fluid flow problem that governs the acoustic phenomena. 
If P* is the ambient pressure and h is the scale height, then the nondimensional form of the Euler 
equations (the equations of continuity, balance of momentum, and energy) is 


in 

at 


+ div(pq) = 0 


+ (q • V )q = — — V p - -Ik 

at yp y 

* N 

i +q V 


(A.l) 
(A. 2) 

(A. 3) 


_ i „ 

Note that in equation (A. 2) the forcing term -^-k 


(k is the unit vector in the z direction) arises due 


to the forcing term per unit mass -gk in the original variables which is due to gravity. In equation 
(A. 3), f(x,y,z,t) dictates the space time dependency of the source, and e measures the energy release 
per unit volume. For the case of an instantaneous energy release, e is given by 


e = 


(Y ~ l)Qo 
h 3 P* 


where Q 0 is the total energy released at time t = 0. The initial conditions are 


(A. 4) 


p=p=e z , q = 0 at t = 0, (A. 5) 

which represent a calm atmosphere and exponentially decaying pressure and density. 

The boundary condition at z = 0 is 


qz = 0, (A. 6) 

which states that the vertical component of the flow is zero at z = 0. 

The acoustic expansion is based on e << 1 and represents the flow as small changes superim- 
posed on the flow of the ambient state. We note that the ambient velocity is zero, but pressure and 
density have the form e _z . Thus, the expansions are 
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q = eu + eV + • • • (A.7) 

p = e _z [l + ep + e 2 p! + • • • j (A. 8) 

p = e -z ^1 + ea + e 2 a x + • • • j (A. 9) 

where u = (u,w) and Uj = (u^wO. Quantities u, U! and w, Wj are the r and z components of the 
acoustic velocities, respectively in cylindrical coordinates. We substitute expansions (A.7)-(A.9) into 
equations (A.1)-(A.3), initial conditions (A. 5), and boundary conditions (A. 6) and retain terms of 
order e to obtain the field equations. 


Time History of P 
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5 















- 24 - 


References 

[1] Cole, J.D. and C. Greifinger, "Acoustic Gravity Waves Produced by Energy Release," RAND 
Corp. Report RM-5738, 1968. 

[2] Hariharan, S.I., "Nonlinear Acoustic Wave Propagation in Atmosphere," ICASE Report No. 86- 
10, 1986. 

[3] Hariharan, S.I., "A Model Problem for Acoustic Wave Propagation in the Atmosphere.” Proceed- 
ings of the First IMACS Symposium on Computational Acoustics, North Holland, Eds. D. Lee, R.L. 
Sternberg, and M.H. Schultz (in press). 

[4] Hariharan, S.I. and H.C. Lester, "Acoustic Shocks in a Variable Area Duct Containing Near Sonic 
Flows,"/. Comp. Phys., Vol. 58, No. 1 (1985). 

[5] Engquist, B. and A. Majda, "Absorbing Boundary Conditions for Acoustic and Elastic Wave Cal- 
culations," Comm. Pure Appl. Math., Vol. 32, pp. 312-358, 1979. 

[6] Bamberger, A., L. Halpern, P. Jolly, and B. Engquist, "The Paraxial Approximation for the Wave 
Equation: Some New Results," Advances in Computer Methods for Partial Differential Equations - V, 
(IMACS), Eds. R. Vichnevetsky and R.S. Stapleman (1984). 

[7] Kreiss, H.O., "Initial Boundary Value Problems for Hyperbolic Systems," Comm. Pure Appl. Math., 
Vol. 23, pp. 277-298, 1970. 

[8] Abarbanel, S. and D. Gottlieb, "Optimal time splitting for two and three dimensional Navier- 
Stokes equations with mixed derivatives," /. Comp. Phys., Vol. 41, No. 1 (May 1981), pp. 1-33. 


Standard Bibliographic Page 


X. Report No. NASA CR- 178270 2. Government Accession No. 

ICASE Report No. 87-19 

3. Recipient’s Catalog No. 

4. Title and Subtitle 

ACOUSTIC GRAVITY WAVES: A COMPUTATIONAL APPROACH 

5. Report Date 

March 1987 

6. Performing Organization Code 

7. Author(s) 

S. I. Hariharan and P. K. Dutt 

8. Performing Organization Report No. 

87-19 

9. Applications in Science 

and Engineering 

Mail Stop 132C, NASA Langley Research Center 
Hampton, VA 23665-5225 

10. Work Unit No. 

“• c»i^?5T n° 

13. Type of Report and Period Covered 

Contractor Report 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 

14. Sponsoring Agency Code 

505-90-21-01 


15. Supplementary Notes 


Langley Technical Monitor: Submitted to Applied Numerical 

J. C. South Mathematics 

Final Report 
16. Abstract 

This paper discusses numerical solutions of a hyperbolic initial boundary 
value problem that arises from acoustic wave propagation in the atmosphere. 
Field equations are derived from the atmospheric fluid flow governed by the 
Euler equations. The resulting original problem is nonlinear. A first order 
linearized version of the problem is used for computational purposes. The main 
difficulty in the problem as with any open boundary problem is in obtaining 
stable boundary conditions. Approximate boundary conditions are derived and 
shown to be stable. Numerical results are presented to verify the effectiveness 
of these boundary conditions. 


17. Key Words (Suggested by Authors(s)) 

atmospheric acoustics, hyperbolic 

1 18. Distribution Statement 

64 - Numerical Analysis 

systems, well-posedness , finite 
differences 


71 - Acoustics 
Unclassified - unlimited 

19. Security Classif.fof this report) 
Unclassified 

20. ^jfecuyty this page) 

21. of Pages 

22 ‘ Ad? 


For sale by the National Technical Information Service, Springfield, Virginia 22161 


NASA-Langley, 1987 




